/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of single_lambda_emergence methods.

#include "single_lambda_emergence.h"
#include "emergence-fits.h"
#include "CCfits/CCfits"
#include "sstream"
#include "biniostream.h"
#include "fits-utilities.h"
#include "constants.h"

using namespace CCfits;
using namespace blitz;

/** Constructor sets up the emergence object from a FITS file.  All
    functionality is in the emergence class.  */
mcrx::single_lambda_emergence::single_lambda_emergence (const CCfits::FITS& file):
    emergence<T_float> (file)
{};

/** Creates FITS data cubes for saving the camera image arrays.  The
    number of slices (i.e., wavelengths) is specified. \bug Doesn't
    work for 1 slice. 

    Also converts the image units. The image units are per area in the
    length units of the camera distance units, which needs to be
    converted to m. We also need to set the solid angle
    normalization. The camera response function was used when adding
    rays. This response function is dA/dOmega for the position on the
    image plane, where A is the area on the normalized image
    plane. What remains is to divide out the area of one pixel on the
    normalized image plane, which is just 1/(nx*ny).  */
void mcrx::single_lambda_emergence::create_HDU (FITS& file, int slices,
						const T_unit_map& units) const
{
  // convert image values (which is flux in the pixels with length
  // units the same as the camera distance units) to surface
  // brightness. This involves converting the 1/area unit to 1/m and
  // dividing by the solid angle of the pixels. We use units to
  // convert whatever length unit we have to m.
  const T_float to_m = units::convert(units.get("length"),"m");
  const string sb_unit = units.get("L_lambda") + "/m^2/sr";
  const string sb_factr_unit = "m^2/(" + units.get("length") + "^2*sr)";
  // (1kpc^2 = 9.5214e38 m^2)

  for (int c=0; c< n_cameras(); ++c) {
    const T_camera& cam = get_camera (c);

    const array_2& image = cam.get_image ();
    assert (image.size()> 0);
    std::ostringstream ost;
    ost << "CAMERA" << c << suffix (); 

    // We only want to create the HDU's if they don't already exist,
    ExtHDU* hdu;
    try {
      hdu = &open_HDU (file,ost.str());
    }
    catch (CCfits::FITS::NoSuchHDU&) {
      std::vector<long> naxes;
      naxes.push_back(image.extent(firstDim ) );
      naxes.push_back(image.extent(secondDim ) );
      naxes.push_back(slices);
      hdu = file.addImage (ost.str(), FLOAT_IMG, naxes);
      //std::cout << "creating HDU " << ost.str()<< endl;

      hdu->writeComment(comment ());

      // write wcs keywords
      cam.write_wcs(*hdu, units);

      const T_float sb_factr = 1./(to_m*to_m*cam.pixel_normalized_area());
      assert(sb_factr > 0);

      // write units
      hdu->addKey("imunit", sb_unit,
		  "Pixel quantity is surface brightness");
      hdu->addKey("sb_factr", sb_factr,"["+sb_factr_unit +
		  "] Internal camera flux to surface brightness");
    }
  } 
  
}

/** Writes the images to a slice in the data cubes in HDU's
    CAMERAi-suffix ().  The HDU's must have been created with a
    previous call to write_HDU. For the file format specification see
    the mcrx-readme. */
void mcrx::single_lambda_emergence::write_images (CCfits::FITS& file,
						  T_float normalization,
						  int slice) const
{
  assert(normalization > 0);

  cout << "Writing " << n_cameras() << " cameras: \n";

  for (int c=0; c< n_cameras(); ++c) {
    // find HDU
    std::ostringstream ost;
    ost << "CAMERA" << c << suffix ();
    cout << ost.str() << endl;

    ExtHDU& hdu = open_HDU (file, ost.str ());

    std::ostringstream ost2;
    ost2 << "normalization" << slice;
    hdu.addKey(ost2.str(), normalization,
	       "Image normalization (internal usage)");

    double sb_factr = 1; // for some reason CCfits makes a double keyword!?
    try {
      hdu.readKey ("sb_factr", sb_factr);
    }
    catch (CCfits::HDU::NoSuchKeyword&) {}

    // Write image array 
    // We have storage order issues here, FITS images are column major
    const array_2& image = get_camera (c).get_image();
    assert (all (image < blitz::huge (T_float ())));
    assert (all (image >= 0));
    array_2 temp_image (image.shape(), ColumnMajorArray<2> ());
    assert(temp_image.isStorageContiguous());

    temp_image = (normalization*sb_factr)*image;
    std::valarray<T_float> temp (temp_image.dataFirst(), temp_image.size());

    // write the proper slice of the data cube
    std::vector<long> begin, end;
    begin.push_back(1 );
    begin.push_back(1 );
    begin.push_back(slice + 1);
    end.push_back(temp_image.extent(firstDim ) );
    end.push_back(temp_image.extent(secondDim) );
    end.push_back(slice + 1  );
    int s = 1;
    for (vector<long>::const_iterator ii = begin.begin(), jj = end.begin(); 
	 ii != begin.end(); ++ii, ++jj) {
      s*= (*jj - *ii + 1);
    }
    assert (s == temp.size());

    hdu.write(begin, end, temp );
  }
}

/** Loads the camera image arrays from a previously saved FITS file. */    
void mcrx::single_lambda_emergence::load_images (const CCfits::FITS& file,
						 int slice)
{
  cout << "Loading " << n_cameras() << " cameras: \n";
  for (int c=0; c< n_cameras(); ++c) {
    T_camera& cam = get_camera (c);

    std::ostringstream ost;
    ost << "CAMERA" << c << suffix ();
    cout << ost.str() << endl;
    //const_cast<CCfits::FITS&> (file).read(string (ost.str()));
    const CCfits::ExtHDU& hdu = open_HDU (file,ost.str() );
    assert (hdu.axis(0) == cam.xsize ());
    assert (hdu.axis(1 ) == cam.ysize ());
    array_2& image = cam.get_image ();

    //read the proper slice of the data cube
    std::vector<long> begin, end;
    begin.push_back(1 );
    begin.push_back(1 );
    begin.push_back(slice + 1);
    end.push_back(image.extent(firstDim ) );
    end.push_back(image.extent(secondDim) );
    end.push_back(slice + 1  );
    int s = 1;
    for (vector<long>::const_iterator ii = begin.begin(), jj = end.begin(); 
	 ii != begin.end(); ++ii, ++jj) {
      s*= (*jj -*ii + 1);
    }

    read(hdu, image);

    // check for a normalization keyword
    float normalization = 1 ;
    double sb_factr = 1;
    try {
      const_cast<CCfits::ExtHDU&> (hdu).readKey (string ("normalization"),
						 normalization);
    }
    catch (...) {}// CCfits::HDU::NoSuchKeyword&) {}
    try {
      const_cast<CCfits::ExtHDU&> (hdu).readKey (string ("sb_factr"),
						sb_factr);
    }
    catch (...) {}// CCfits::HDU::NoSuchKeyword&) {}

    image *= 1/(normalization*sb_factr);
  }
}


/** Dump the image arrays to a file.  The main objective here is to be
    fast, since we have gotten a termination signal. */
void mcrx::single_lambda_emergence::write_dump (binofstream& file) const
{
  for (int c=0; c< n_cameras(); ++c) {
    const array_2& image = get_camera (c).get_image();
    // We strictly don't need to save the image size here but it's a good check
    array_2::T_index extent = image.shape();
    file << extent [0]
	 << extent [1]; 

    assert(image.isStorageContiguous());
    file.write(reinterpret_cast<const char*> (image.dataFirst()),
	       image.storageSize()*sizeof (array_2::T_numtype) );
  } 
}

/** Read the image arrays from a previously written dump file. */
void mcrx::single_lambda_emergence::load_dump (binifstream& file)
{
  for (int c=0; c< n_cameras(); ++c) {
    T_camera& cam = get_camera (c);

    array_2& image = cam.get_image ();
    array_2::T_index size;
    file >> size [0]
	 >>  size [1]; 
    //cout << size << endl;
    assert (size [0] == cam.xsize ());
    assert (size [1] == cam.ysize ());

    assert(image.isStorageContiguous());
    file.read (reinterpret_cast<char*> (image.dataFirst()),
	       image.storageSize()*sizeof (array_2::T_numtype) );
  }
}


